!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!

! **************************************************************************************************
!> \brief Working with the DFTB parameter types.
!> \author JGH (24.02.2007)
! **************************************************************************************************
MODULE qs_dftb_utils

   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
                                              cp_logger_type
   USE cp_output_handling,              ONLY: cp_p_file,&
                                              cp_print_key_finished_output,&
                                              cp_print_key_should_output,&
                                              cp_print_key_unit_nr
   USE input_section_types,             ONLY: section_vals_type
   USE kinds,                           ONLY: default_string_length,&
                                              dp
   USE qs_dftb_types,                   ONLY: qs_dftb_atom_type
#include "./base/base_uses.f90"

   IMPLICIT NONE

   PRIVATE

   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dftb_utils'

   ! Maximum number of points used for interpolation
   INTEGER, PARAMETER                     :: max_inter = 5
   ! Maximum number of points used for extrapolation
   INTEGER, PARAMETER                     :: max_extra = 9
   ! see also qs_dftb_parameters
   REAL(dp), PARAMETER                    :: slako_d0 = 1._dp
   ! pointer to skab
   INTEGER, DIMENSION(0:3, 0:3, 0:3, 0:3, 0:3):: iptr
   ! small real number
   REAL(dp), PARAMETER                    :: rtiny = 1.e-10_dp
   ! eta(0) for mm atoms and non-scc qm atoms
   REAL(dp), PARAMETER                    :: eta_mm = 0.47_dp
   ! step size for qmmm finite difference
   REAL(dp), PARAMETER                    :: ddrmm = 0.0001_dp

   PUBLIC :: allocate_dftb_atom_param, &
             deallocate_dftb_atom_param, &
             get_dftb_atom_param, &
             set_dftb_atom_param, &
             write_dftb_atom_param
   PUBLIC :: compute_block_sk, &
             urep_egr, iptr

CONTAINS

! **************************************************************************************************
!> \brief ...
!> \param dftb_parameter ...
! **************************************************************************************************
   SUBROUTINE allocate_dftb_atom_param(dftb_parameter)

      TYPE(qs_dftb_atom_type), POINTER                   :: dftb_parameter

      IF (ASSOCIATED(dftb_parameter)) &
         CALL deallocate_dftb_atom_param(dftb_parameter)

      ALLOCATE (dftb_parameter)

      dftb_parameter%defined = .FALSE.
      dftb_parameter%name = ""
      dftb_parameter%typ = "NONE"
      dftb_parameter%z = -1
      dftb_parameter%zeff = -1.0_dp
      dftb_parameter%natorb = 0
      dftb_parameter%lmax = -1
      dftb_parameter%skself = 0.0_dp
      dftb_parameter%occupation = 0.0_dp
      dftb_parameter%eta = 0.0_dp
      dftb_parameter%energy = 0.0_dp
      dftb_parameter%xi = 0.0_dp
      dftb_parameter%di = 0.0_dp
      dftb_parameter%rcdisp = 0.0_dp
      dftb_parameter%dudq = 0.0_dp

   END SUBROUTINE allocate_dftb_atom_param

! **************************************************************************************************
!> \brief ...
!> \param dftb_parameter ...
! **************************************************************************************************
   SUBROUTINE deallocate_dftb_atom_param(dftb_parameter)

      TYPE(qs_dftb_atom_type), POINTER                   :: dftb_parameter

      CPASSERT(ASSOCIATED(dftb_parameter))
      DEALLOCATE (dftb_parameter)

   END SUBROUTINE deallocate_dftb_atom_param

! **************************************************************************************************
!> \brief ...
!> \param dftb_parameter ...
!> \param name ...
!> \param typ ...
!> \param defined ...
!> \param z ...
!> \param zeff ...
!> \param natorb ...
!> \param lmax ...
!> \param skself ...
!> \param occupation ...
!> \param eta ...
!> \param energy ...
!> \param cutoff ...
!> \param xi ...
!> \param di ...
!> \param rcdisp ...
!> \param dudq ...
! **************************************************************************************************
   SUBROUTINE get_dftb_atom_param(dftb_parameter, name, typ, defined, z, zeff, natorb, &
                                  lmax, skself, occupation, eta, energy, cutoff, xi, di, rcdisp, dudq)

      TYPE(qs_dftb_atom_type), POINTER                   :: dftb_parameter
      CHARACTER(LEN=default_string_length), &
         INTENT(OUT), OPTIONAL                           :: name, typ
      LOGICAL, INTENT(OUT), OPTIONAL                     :: defined
      INTEGER, INTENT(OUT), OPTIONAL                     :: z
      REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: zeff
      INTEGER, INTENT(OUT), OPTIONAL                     :: natorb, lmax
      REAL(KIND=dp), DIMENSION(0:3), OPTIONAL            :: skself, occupation, eta
      REAL(KIND=dp), OPTIONAL                            :: energy, cutoff, xi, di, rcdisp, dudq

      CPASSERT(ASSOCIATED(dftb_parameter))

      IF (PRESENT(name)) name = dftb_parameter%name
      IF (PRESENT(typ)) typ = dftb_parameter%typ
      IF (PRESENT(defined)) defined = dftb_parameter%defined
      IF (PRESENT(z)) z = dftb_parameter%z
      IF (PRESENT(zeff)) zeff = dftb_parameter%zeff
      IF (PRESENT(natorb)) natorb = dftb_parameter%natorb
      IF (PRESENT(lmax)) lmax = dftb_parameter%lmax
      IF (PRESENT(skself)) skself = dftb_parameter%skself
      IF (PRESENT(eta)) eta = dftb_parameter%eta
      IF (PRESENT(energy)) energy = dftb_parameter%energy
      IF (PRESENT(cutoff)) cutoff = dftb_parameter%cutoff
      IF (PRESENT(occupation)) occupation = dftb_parameter%occupation
      IF (PRESENT(xi)) xi = dftb_parameter%xi
      IF (PRESENT(di)) di = dftb_parameter%di
      IF (PRESENT(rcdisp)) rcdisp = dftb_parameter%rcdisp
      IF (PRESENT(dudq)) dudq = dftb_parameter%dudq

   END SUBROUTINE get_dftb_atom_param

! **************************************************************************************************
!> \brief ...
!> \param dftb_parameter ...
!> \param name ...
!> \param typ ...
!> \param defined ...
!> \param z ...
!> \param zeff ...
!> \param natorb ...
!> \param lmax ...
!> \param skself ...
!> \param occupation ...
!> \param eta ...
!> \param energy ...
!> \param cutoff ...
!> \param xi ...
!> \param di ...
!> \param rcdisp ...
!> \param dudq ...
! **************************************************************************************************
   SUBROUTINE set_dftb_atom_param(dftb_parameter, name, typ, defined, z, zeff, natorb, &
                                  lmax, skself, occupation, eta, energy, cutoff, xi, di, rcdisp, dudq)

      TYPE(qs_dftb_atom_type), POINTER                   :: dftb_parameter
      CHARACTER(LEN=default_string_length), INTENT(IN), &
         OPTIONAL                                        :: name, typ
      LOGICAL, INTENT(IN), OPTIONAL                      :: defined
      INTEGER, INTENT(IN), OPTIONAL                      :: z
      REAL(KIND=dp), INTENT(IN), OPTIONAL                :: zeff
      INTEGER, INTENT(IN), OPTIONAL                      :: natorb, lmax
      REAL(KIND=dp), DIMENSION(0:3), OPTIONAL            :: skself, occupation, eta
      REAL(KIND=dp), OPTIONAL                            :: energy, cutoff, xi, di, rcdisp, dudq

      CPASSERT(ASSOCIATED(dftb_parameter))

      IF (PRESENT(name)) dftb_parameter%name = name
      IF (PRESENT(typ)) dftb_parameter%typ = typ
      IF (PRESENT(defined)) dftb_parameter%defined = defined
      IF (PRESENT(z)) dftb_parameter%z = z
      IF (PRESENT(zeff)) dftb_parameter%zeff = zeff
      IF (PRESENT(natorb)) dftb_parameter%natorb = natorb
      IF (PRESENT(lmax)) dftb_parameter%lmax = lmax
      IF (PRESENT(skself)) dftb_parameter%skself = skself
      IF (PRESENT(eta)) dftb_parameter%eta = eta
      IF (PRESENT(occupation)) dftb_parameter%occupation = occupation
      IF (PRESENT(energy)) dftb_parameter%energy = energy
      IF (PRESENT(cutoff)) dftb_parameter%cutoff = cutoff
      IF (PRESENT(xi)) dftb_parameter%xi = xi
      IF (PRESENT(di)) dftb_parameter%di = di
      IF (PRESENT(rcdisp)) dftb_parameter%rcdisp = rcdisp
      IF (PRESENT(dudq)) dftb_parameter%dudq = dudq

   END SUBROUTINE set_dftb_atom_param

! **************************************************************************************************
!> \brief ...
!> \param dftb_parameter ...
!> \param subsys_section ...
! **************************************************************************************************
   SUBROUTINE write_dftb_atom_param(dftb_parameter, subsys_section)

      TYPE(qs_dftb_atom_type), POINTER                   :: dftb_parameter
      TYPE(section_vals_type), POINTER                   :: subsys_section

      CHARACTER(LEN=default_string_length)               :: name, typ
      INTEGER                                            :: lmax, natorb, output_unit, z
      LOGICAL                                            :: defined
      REAL(dp)                                           :: zeff
      TYPE(cp_logger_type), POINTER                      :: logger

      NULLIFY (logger)
      logger => cp_get_default_logger()
      IF (ASSOCIATED(dftb_parameter) .AND. &
          BTEST(cp_print_key_should_output(logger%iter_info, subsys_section, &
                                           "PRINT%KINDS/POTENTIAL"), cp_p_file)) THEN

         output_unit = cp_print_key_unit_nr(logger, subsys_section, "PRINT%KINDS", &
                                            extension=".Log")

         IF (output_unit > 0) THEN
            CALL get_dftb_atom_param(dftb_parameter, name=name, typ=typ, defined=defined, &
                                     z=z, zeff=zeff, natorb=natorb, lmax=lmax)

            WRITE (UNIT=output_unit, FMT="(/,A,T67,A14)") &
               " DFTB  parameters: ", TRIM(name)
            IF (defined) THEN
               WRITE (UNIT=output_unit, FMT="(T16,A,T71,F10.2)") &
                  "Effective core charge:", zeff
               WRITE (UNIT=output_unit, FMT="(T16,A,T71,I10)") &
                  "Number of orbitals:", natorb
            ELSE
               WRITE (UNIT=output_unit, FMT="(T55,A)") &
                  "Parameters are not defined"
            END IF
         END IF
         CALL cp_print_key_finished_output(output_unit, logger, subsys_section, &
                                           "PRINT%KINDS")
      END IF

   END SUBROUTINE write_dftb_atom_param

! **************************************************************************************************
!> \brief ...
!> \param block ...
!> \param smatij ...
!> \param smatji ...
!> \param rij ...
!> \param ngrd ...
!> \param ngrdcut ...
!> \param dgrd ...
!> \param llm ...
!> \param lmaxi ...
!> \param lmaxj ...
!> \param irow ...
!> \param iatom ...
! **************************************************************************************************
   SUBROUTINE compute_block_sk(block, smatij, smatji, rij, ngrd, ngrdcut, dgrd, &
                               llm, lmaxi, lmaxj, irow, iatom)
      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: block, smatij, smatji
      REAL(KIND=dp), DIMENSION(3)                        :: rij
      INTEGER                                            :: ngrd, ngrdcut
      REAL(KIND=dp)                                      :: dgrd
      INTEGER                                            :: llm, lmaxi, lmaxj, irow, iatom

      REAL(KIND=dp)                                      :: dr
      REAL(KIND=dp), DIMENSION(20)                       :: skabij, skabji

      dr = SQRT(SUM(rij(:)**2))
      CALL getskz(smatij, skabij, dr, ngrd, ngrdcut, dgrd, llm)
      CALL getskz(smatji, skabji, dr, ngrd, ngrdcut, dgrd, llm)
      IF (irow == iatom) THEN
         CALL turnsk(block, skabji, skabij, rij, dr, lmaxi, lmaxj)
      ELSE
         CALL turnsk(block, skabij, skabji, -rij, dr, lmaxj, lmaxi)
      END IF

   END SUBROUTINE compute_block_sk

! **************************************************************************************************
!> \brief Gets matrix elements on z axis, as they are stored in the tables
!> \param slakotab ...
!> \param skpar ...
!> \param dx ...
!> \param ngrd ...
!> \param ngrdcut ...
!> \param dgrd ...
!> \param llm ...
!> \author 07. Feb. 2004, TH
! **************************************************************************************************
   SUBROUTINE getskz(slakotab, skpar, dx, ngrd, ngrdcut, dgrd, llm)
      REAL(dp), INTENT(in)                               :: slakotab(:, :), dx
      INTEGER, INTENT(in)                                :: ngrd, ngrdcut
      REAL(dp), INTENT(in)                               :: dgrd
      INTEGER, INTENT(in)                                :: llm
      REAL(dp), INTENT(out)                              :: skpar(llm)

      INTEGER                                            :: clgp

      skpar = 0._dp
      !
      ! Determine closest grid point
      !
      clgp = NINT(dx/dgrd)
      !
      ! Screen elements which are too far away
      !
      IF (clgp > ngrdcut) RETURN
      !
      ! The grid point is either contained in the table --> matrix element
      ! can be interpolated, or it is outside the table --> matrix element
      ! needs to be extrapolated.
      !
      IF (clgp > ngrd) THEN
         !
         ! Extrapolate external matrix elements if table does not finish with zero
         !
         CALL extrapol(slakotab, skpar, dx, ngrd, dgrd, llm)
      ELSE
         !
         ! Interpolate tabulated matrix elements
         !
         CALL interpol(slakotab, skpar, dx, ngrd, dgrd, llm, clgp)
      END IF
   END SUBROUTINE getskz

! **************************************************************************************************
!> \brief ...
!> \param slakotab ...
!> \param skpar ...
!> \param dx ...
!> \param ngrd ...
!> \param dgrd ...
!> \param llm ...
!> \param clgp ...
! **************************************************************************************************
   SUBROUTINE interpol(slakotab, skpar, dx, ngrd, dgrd, llm, clgp)
      REAL(dp), INTENT(in)                               :: slakotab(:, :), dx
      INTEGER, INTENT(in)                                :: ngrd
      REAL(dp), INTENT(in)                               :: dgrd
      INTEGER, INTENT(in)                                :: llm
      REAL(dp), INTENT(out)                              :: skpar(llm)
      INTEGER, INTENT(in)                                :: clgp

      INTEGER                                            :: fgpm, k, l, lgpm
      REAL(dp)                                           :: error, xa(max_inter), ya(max_inter)

      lgpm = MIN(clgp + INT(max_inter/2.0), ngrd)
      fgpm = lgpm - max_inter + 1
      DO k = 0, max_inter - 1
         xa(k + 1) = (fgpm + k)*dgrd
      END DO
      !
      ! Interpolate matrix elements for all orbitals
      !
      DO l = 1, llm
         !
         ! Read SK parameters from table
         !
         ya(1:max_inter) = slakotab(fgpm:lgpm, l)
         CALL polint(xa, ya, max_inter, dx, skpar(l), error)
      END DO
   END SUBROUTINE interpol

! **************************************************************************************************
!> \brief ...
!> \param slakotab ...
!> \param skpar ...
!> \param dx ...
!> \param ngrd ...
!> \param dgrd ...
!> \param llm ...
! **************************************************************************************************
   SUBROUTINE extrapol(slakotab, skpar, dx, ngrd, dgrd, llm)
      REAL(dp), INTENT(in)                               :: slakotab(:, :), dx
      INTEGER, INTENT(in)                                :: ngrd
      REAL(dp), INTENT(in)                               :: dgrd
      INTEGER, INTENT(in)                                :: llm
      REAL(dp), INTENT(out)                              :: skpar(llm)

      INTEGER                                            :: fgp, k, l, lgp, ntable, nzero
      REAL(dp)                                           :: error, xa(max_extra), ya(max_extra)

      nzero = max_extra/3
      ntable = max_extra - nzero
      !
      ! Get the three last distances from the table
      !
      DO k = 1, ntable
         xa(k) = (ngrd - (max_extra - 3) + k)*dgrd
      END DO
      DO k = 1, nzero
         xa(ntable + k) = (ngrd + k - 1)*dgrd + slako_d0
         ya(ntable + k) = 0.0
      END DO
      !
      ! Extrapolate matrix elements for all orbitals
      !
      DO l = 1, llm
         !
         ! Read SK parameters from table
         !
         fgp = ngrd + 1 - (max_extra - 3)
         lgp = ngrd
         ya(1:max_extra - 3) = slakotab(fgp:lgp, l)
         CALL polint(xa, ya, max_extra, dx, skpar(l), error)
      END DO
   END SUBROUTINE extrapol

! **************************************************************************************************
!> \brief   Turn matrix element from z-axis to orientation of dxv
!> \param mat ...
!> \param skab1 ...
!> \param skab2 ...
!> \param dxv ...
!> \param dx ...
!> \param lmaxa ...
!> \param lmaxb ...
!> \date    13. Jan 2004
!> \par Notes
!>          These routines are taken from an old TB code (unknown to TH).
!>          They are highly optimised and taken because they are time critical.
!>          They are explicit, so not recursive, and work up to d functions.
!>
!>          Set variables necessary for rotation of matrix elements
!>
!>          r_i^2/r, replicated in rr2(4:6) for index convenience later
!>          r_i/r, direction vector, rr(4:6) are replicated from 1:3
!>          lmax of A and B
!> \author  TH
!> \version 1.0
! **************************************************************************************************
   SUBROUTINE turnsk(mat, skab1, skab2, dxv, dx, lmaxa, lmaxb)
      REAL(dp), INTENT(inout)                  :: mat(:, :)
      REAL(dp), INTENT(in)                     :: skab1(:), skab2(:), dxv(3), dx
      INTEGER, INTENT(in)                      :: lmaxa, lmaxb

      INTEGER                                  :: lmaxab, minlmaxab
      REAL(dp)                                 :: rinv, rr(6), rr2(6)

      lmaxab = MAX(lmaxa, lmaxb)
      ! Determine l quantum limits.
      IF (lmaxab .GT. 2) CPABORT('lmax=2')
      minlmaxab = MIN(lmaxa, lmaxb)
      !
      ! s-s interaction
      !
      CALL skss(skab1, mat)
      !
      IF (lmaxab .LE. 0) RETURN
      !
      rr2(1:3) = dxv(1:3)**2
      rr(1:3) = dxv(1:3)
      rinv = 1.0_dp/dx
      !
      rr(1:3) = rr(1:3)*rinv
      rr(4:6) = rr(1:3)
      rr2(1:3) = rr2(1:3)*rinv**2
      rr2(4:6) = rr2(1:3)
      !
      ! s-p, p-s and p-p interaction
      !
      IF (minlmaxab .GE. 1) THEN
         CALL skpp(skab1, mat, iptr(:, :, :, lmaxa, lmaxb))
         CALL sksp(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
         CALL sksp(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.)
      ELSE
         IF (lmaxb .GE. 1) THEN
            CALL sksp(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
         ELSE
            CALL sksp(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.)
         END IF
      END IF
      !
      ! If there is only s-p interaction we have finished
      !
      IF (lmaxab .LE. 1) RETURN
      !
      ! at least one atom has d functions
      !
      IF (minlmaxab .EQ. 2) THEN
         !
         ! in case both atoms have d functions
         !
         CALL skdd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb))
         CALL sksd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
         CALL sksd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.)
         CALL skpd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
         CALL skpd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.)
      ELSE
         !
         ! One atom has d functions, the other has s or s and p functions
         !
         IF (lmaxa .EQ. 0) THEN
            !
            ! atom b has d, the atom a only s functions
            !
            CALL sksd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
         ELSE IF (lmaxa .EQ. 1) THEN
            !
            ! atom b has d, the atom a s and p functions
            !
            CALL sksd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
            CALL skpd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
         ELSE
            !
            ! atom a has d functions
            !
            IF (lmaxb .EQ. 0) THEN
               !
               ! atom a has d, atom b has only s functions
               !
               CALL sksd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.)
            ELSE
               !
               ! atom a has d, atom b has s and p functions
               !
               CALL sksd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.)
               CALL skpd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.)
            END IF
         END IF
      END IF
      !
   CONTAINS
      !
      ! The subroutines to turn the matrix elements are taken as internal subroutines
      ! as it is beneficial to inline them.
      !
      ! They are both turning the matrix elements and placing them appropriately
      ! into the matrix block
      !
! **************************************************************************************************
!> \brief   s-s interaction (no rotation necessary)
!> \param skpar ...
!> \param mat ...
!> \version 1.0
! **************************************************************************************************
      SUBROUTINE skss(skpar, mat)
      REAL(dp), INTENT(in)                               :: skpar(:)
      REAL(dp), INTENT(inout)                            :: mat(:, :)

         mat(1, 1) = mat(1, 1) + skpar(1)
         !
      END SUBROUTINE skss

! **************************************************************************************************
!> \brief  s-p interaction (simple rotation)
!> \param skpar ...
!> \param mat ...
!> \param ind ...
!> \param transposed ...
!> \version 1.0
! **************************************************************************************************
      SUBROUTINE sksp(skpar, mat, ind, transposed)
      REAL(dp), INTENT(in)                               :: skpar(:)
      REAL(dp), INTENT(inout)                            :: mat(:, :)
      INTEGER, INTENT(in)                                :: ind(0:, 0:, 0:)
      LOGICAL, INTENT(in)                                :: transposed

      INTEGER                                            :: l
      REAL(dp)                                           :: skp

         skp = skpar(ind(1, 0, 0))
         IF (transposed) THEN
            DO l = 1, 3
               mat(1, l + 1) = mat(1, l + 1) + rr(l)*skp
            END DO
         ELSE
            DO l = 1, 3
               mat(l + 1, 1) = mat(l + 1, 1) - rr(l)*skp
            END DO
         END IF
         !
      END SUBROUTINE sksp

! **************************************************************************************************
!> \brief ...
!> \param skpar ...
!> \param mat ...
!> \param ind ...
! **************************************************************************************************
      SUBROUTINE skpp(skpar, mat, ind)
      REAL(dp), INTENT(in)                               :: skpar(:)
      REAL(dp), INTENT(inout)                            :: mat(:, :)
      INTEGER, INTENT(in)                                :: ind(0:, 0:, 0:)

      INTEGER                                            :: ii, ir, is, k, l
      REAL(dp)                                           :: epp(6), matel(6), skppp, skpps

         epp(1:3) = rr2(1:3)
         DO l = 1, 3
            epp(l + 3) = rr(l)*rr(l + 1)
         END DO
         skppp = skpar(ind(1, 1, 1))
         skpps = skpar(ind(1, 1, 0))
         !
         DO l = 1, 3
            matel(l) = epp(l)*skpps + (1._dp - epp(l))*skppp
         END DO
         DO l = 4, 6
            matel(l) = epp(l)*(skpps - skppp)
         END DO
         !
         DO ir = 1, 3
            DO is = 1, ir - 1
               ii = ir - is
               k = 3*ii - (ii*(ii - 1))/2 + is
               mat(is + 1, ir + 1) = mat(is + 1, ir + 1) + matel(k)
               mat(ir + 1, is + 1) = mat(ir + 1, is + 1) + matel(k)
            END DO
            mat(ir + 1, ir + 1) = mat(ir + 1, ir + 1) + matel(ir)
         END DO
      END SUBROUTINE skpp

! **************************************************************************************************
!> \brief ...
!> \param skpar ...
!> \param mat ...
!> \param ind ...
!> \param transposed ...
! **************************************************************************************************
      SUBROUTINE sksd(skpar, mat, ind, transposed)
      REAL(dp), INTENT(in)                               :: skpar(:)
      REAL(dp), INTENT(inout)                            :: mat(:, :)
      INTEGER, INTENT(in)                                :: ind(0:, 0:, 0:)
      LOGICAL, INTENT(in)                                :: transposed

      INTEGER                                            :: l
      REAL(dp)                                           :: d4, d5, es(5), r3, sksds

         sksds = skpar(ind(2, 0, 0))
         r3 = SQRT(3._dp)
         d4 = rr2(3) - 0.5_dp*(rr2(1) + rr2(2))
         d5 = rr2(1) - rr2(2)
         !
         DO l = 1, 3
            es(l) = r3*rr(l)*rr(l + 1)
         END DO
         es(4) = 0.5_dp*r3*d5
         es(5) = d4
         !
         IF (transposed) THEN
            DO l = 1, 5
               mat(1, l + 4) = mat(1, l + 4) + es(l)*sksds
            END DO
         ELSE
            DO l = 1, 5
               mat(l + 4, 1) = mat(l + 4, 1) + es(l)*sksds
            END DO
         END IF
      END SUBROUTINE sksd

! **************************************************************************************************
!> \brief ...
!> \param skpar ...
!> \param mat ...
!> \param ind ...
!> \param transposed ...
! **************************************************************************************************
      SUBROUTINE skpd(skpar, mat, ind, transposed)
      REAL(dp), INTENT(in)                               :: skpar(:)
      REAL(dp), INTENT(inout)                            :: mat(:, :)
      INTEGER, INTENT(in)                                :: ind(0:, 0:, 0:)
      LOGICAL, INTENT(in)                                :: transposed

      INTEGER                                            :: ir, is, k, l, m
      REAL(dp)                                           :: d3, d4, d5, d6, dm(15), epd(13, 2), r3, &
                                                            sktmp

         r3 = SQRT(3.0_dp)
         d3 = rr2(1) + rr2(2)
         d4 = rr2(3) - 0.5_dp*d3
         d5 = rr2(1) - rr2(2)
         d6 = rr(1)*rr(2)*rr(3)
         DO l = 1, 3
            epd(l, 1) = r3*rr2(l)*rr(l + 1)
            epd(l, 2) = rr(l + 1)*(1.0_dp - 2._dp*rr2(l))
            epd(l + 4, 1) = r3*rr2(l)*rr(l + 2)
            epd(l + 4, 2) = rr(l + 2)*(1.0_dp - 2*rr2(l))
            epd(l + 7, 1) = 0.5_dp*r3*rr(l)*d5
            epd(l + 10, 1) = rr(l)*d4
         END DO
         !
         epd(4, 1) = r3*d6
         epd(4, 2) = -2._dp*d6
         epd(8, 2) = rr(1)*(1.0_dp - d5)
         epd(9, 2) = -rr(2)*(1.0_dp + d5)
         epd(10, 2) = -rr(3)*d5
         epd(11, 2) = -r3*rr(1)*rr2(3)
         epd(12, 2) = -r3*rr(2)*rr2(3)
         epd(13, 2) = r3*rr(3)*d3
         !
         dm(1:15) = 0.0_dp
         !
         DO m = 1, 2
            sktmp = skpar(ind(2, 1, m - 1))
            dm(1) = dm(1) + epd(1, m)*sktmp
            dm(2) = dm(2) + epd(6, m)*sktmp
            dm(3) = dm(3) + epd(4, m)*sktmp
            dm(5) = dm(5) + epd(2, m)*sktmp
            dm(6) = dm(6) + epd(7, m)*sktmp
            dm(7) = dm(7) + epd(5, m)*sktmp
            dm(9) = dm(9) + epd(3, m)*sktmp
            DO l = 8, 13
               dm(l + 2) = dm(l + 2) + epd(l, m)*sktmp
            END DO
         END DO
         !
         dm(4) = dm(3)
         dm(8) = dm(3)
         !
         IF (transposed) THEN
            DO ir = 1, 5
               DO is = 1, 3
                  k = 3*(ir - 1) + is
                  mat(is + 1, ir + 4) = mat(is + 1, ir + 4) + dm(k)
               END DO
            END DO
         ELSE
            DO ir = 1, 5
               DO is = 1, 3
                  k = 3*(ir - 1) + is
                  mat(ir + 4, is + 1) = mat(ir + 4, is + 1) - dm(k)
               END DO
            END DO
         END IF
         !
      END SUBROUTINE skpd

! **************************************************************************************************
!> \brief ...
!> \param skpar ...
!> \param mat ...
!> \param ind ...
! **************************************************************************************************
      SUBROUTINE skdd(skpar, mat, ind)
      REAL(dp), INTENT(in)                               :: skpar(:)
      REAL(dp), INTENT(inout)                            :: mat(:, :)
      INTEGER, INTENT(in)                                :: ind(0:, 0:, 0:)

      INTEGER                                            :: ii, ir, is, k, l, m
      REAL(dp)                                           :: d3, d4, d5, dd(3), dm(15), e(15, 3), r3

         r3 = SQRT(3._dp)
         d3 = rr2(1) + rr2(2)
         d4 = rr2(3) - 0.5_dp*d3
         d5 = rr2(1) - rr2(2)
         DO l = 1, 3
            e(l, 1) = rr2(l)*rr2(l + 1)
            e(l, 2) = rr2(l) + rr2(l + 1) - 4._dp*e(l, 1)
            e(l, 3) = rr2(l + 2) + e(l, 1)
            e(l, 1) = 3._dp*e(l, 1)
         END DO
         e(4, 1) = d5**2
         e(4, 2) = d3 - e(4, 1)
         e(4, 3) = rr2(3) + 0.25_dp*e(4, 1)
         e(4, 1) = 0.75_dp*e(4, 1)
         e(5, 1) = d4**2
         e(5, 2) = 3._dp*rr2(3)*d3
         e(5, 3) = 0.75_dp*d3**2
         dd(1) = rr(1)*rr(3)
         dd(2) = rr(2)*rr(1)
         dd(3) = rr(3)*rr(2)
         DO l = 1, 2
            e(l + 5, 1) = 3._dp*rr2(l + 1)*dd(l)
            e(l + 5, 2) = dd(l)*(1._dp - 4._dp*rr2(l + 1))
            e(l + 5, 3) = dd(l)*(rr2(l + 1) - 1._dp)
         END DO
         e(8, 1) = dd(1)*d5*1.5_dp
         e(8, 2) = dd(1)*(1.0_dp - 2.0_dp*d5)
         e(8, 3) = dd(1)*(0.5_dp*d5 - 1.0_dp)
         e(9, 1) = d5*0.5_dp*d4*r3
         e(9, 2) = -d5*rr2(3)*r3
         e(9, 3) = d5*0.25_dp*(1.0_dp + rr2(3))*r3
         e(10, 1) = rr2(1)*dd(3)*3.0_dp
         e(10, 2) = (0.25_dp - rr2(1))*dd(3)*4.0_dp
         e(10, 3) = dd(3)*(rr2(1) - 1.0_dp)
         e(11, 1) = 1.5_dp*dd(3)*d5
         e(11, 2) = -dd(3)*(1.0_dp + 2.0_dp*d5)
         e(11, 3) = dd(3)*(1.0_dp + 0.5_dp*d5)
         e(13, 3) = 0.5_dp*d5*dd(2)
         e(13, 2) = -2.0_dp*dd(2)*d5
         e(13, 1) = e(13, 3)*3.0_dp
         e(12, 1) = d4*dd(1)*r3
         e(14, 1) = d4*dd(3)*r3
         e(15, 1) = d4*dd(2)*r3
         e(15, 2) = -2.0_dp*r3*dd(2)*rr2(3)
         e(15, 3) = 0.5_dp*r3*(1.0_dp + rr2(3))*dd(2)
         e(14, 2) = r3*dd(3)*(d3 - rr2(3))
         e(14, 3) = -r3*0.5_dp*dd(3)*d3
         e(12, 2) = r3*dd(1)*(d3 - rr2(3))
         e(12, 3) = -r3*0.5_dp*dd(1)*d3
         !
         dm(1:15) = 0._dp
         DO l = 1, 15
            DO m = 1, 3
               dm(l) = dm(l) + e(l, m)*skpar(ind(2, 2, m - 1))
            END DO
         END DO
         !
         DO ir = 1, 5
            DO is = 1, ir - 1
               ii = ir - is
               k = 5*ii - (ii*(ii - 1))/2 + is
               mat(ir + 4, is + 4) = mat(ir + 4, is + 4) + dm(k)
               mat(is + 4, ir + 4) = mat(is + 4, ir + 4) + dm(k)
            END DO
            mat(ir + 4, ir + 4) = mat(ir + 4, ir + 4) + dm(ir)
         END DO
      END SUBROUTINE skdd
      !
   END SUBROUTINE turnsk

! **************************************************************************************************
!> \brief ...
!> \param xa ...
!> \param ya ...
!> \param n ...
!> \param x ...
!> \param y ...
!> \param dy ...
! **************************************************************************************************
   SUBROUTINE polint(xa, ya, n, x, y, dy)
      INTEGER, INTENT(in)                                :: n
      REAL(dp), INTENT(in)                               :: ya(n), xa(n), x
      REAL(dp), INTENT(out)                              :: y, dy

      INTEGER                                            :: i, m, ns
      REAL(dp)                                           :: c(n), d(n), den, dif, dift, ho, hp, w

!
!

      ns = 1

      dif = ABS(x - xa(1))
      DO i = 1, n
         dift = ABS(x - xa(i))
         IF (dift .LT. dif) THEN
            ns = i
            dif = dift
         END IF
         c(i) = ya(i)
         d(i) = ya(i)
      END DO
      !
      y = ya(ns)
      ns = ns - 1
      DO m = 1, n - 1
         DO i = 1, n - m
            ho = xa(i) - x
            hp = xa(i + m) - x
            w = c(i + 1) - d(i)
            den = ho - hp
            CPASSERT(den /= 0.0_dp)
            den = w/den
            d(i) = hp*den
            c(i) = ho*den
         END DO
         IF (2*ns .LT. n - m) THEN
            dy = c(ns + 1)
         ELSE
            dy = d(ns)
            ns = ns - 1
         END IF
         y = y + dy
      END DO
!
      RETURN
   END SUBROUTINE polint

! **************************************************************************************************
!> \brief ...
!> \param rv ...
!> \param r ...
!> \param erep ...
!> \param derep ...
!> \param n_urpoly ...
!> \param urep ...
!> \param spdim ...
!> \param s_cut ...
!> \param srep ...
!> \param spxr ...
!> \param scoeff ...
!> \param surr ...
!> \param dograd ...
! **************************************************************************************************
   SUBROUTINE urep_egr(rv, r, erep, derep, &
                       n_urpoly, urep, spdim, s_cut, srep, spxr, scoeff, surr, dograd)

      REAL(dp), INTENT(in)                               :: rv(3), r
      REAL(dp), INTENT(inout)                            :: erep, derep(3)
      INTEGER, INTENT(in)                                :: n_urpoly
      REAL(dp), INTENT(in)                               :: urep(:)
      INTEGER, INTENT(in)                                :: spdim
      REAL(dp), INTENT(in)                               :: s_cut, srep(3)
      REAL(dp), POINTER                                  :: spxr(:, :), scoeff(:, :)
      REAL(dp), INTENT(in)                               :: surr(2)
      LOGICAL, INTENT(in)                                :: dograd

      INTEGER                                            :: ic, isp, jsp, nsp
      REAL(dp)                                           :: de_z, rz

      derep = 0._dp
      de_z = 0._dp
      IF (n_urpoly > 0) THEN
         !
         ! polynomial part
         !
         rz = urep(1) - r
         IF (rz <= rtiny) RETURN
         DO ic = 2, n_urpoly
            erep = erep + urep(ic)*rz**(ic)
         END DO
         IF (dograd) THEN
            DO ic = 2, n_urpoly
               de_z = de_z - ic*urep(ic)*rz**(ic - 1)
            END DO
         END IF
      ELSE IF (spdim > 0) THEN
         !
         ! spline part
         !
         ! This part is kind of proprietary Paderborn code and I won't reverse-engineer
         ! everything in detail. What is obvious is documented.
         !
         ! This part has 4 regions:
         ! a) very long range is screened
         ! b) short-range is extrapolated with e-functions
         ! ca) normal range is approximated with a spline
         ! cb) longer range is extrapolated with an higher degree spline
         !
         IF (r > s_cut) RETURN ! screening (condition a)
         !
         IF (r < spxr(1, 1)) THEN
            ! a) short range
            erep = erep + EXP(-srep(1)*r + srep(2)) + srep(3)
            IF (dograd) de_z = de_z - srep(1)*EXP(-srep(1)*r + srep(2))
         ELSE
            !
            ! condition c). First determine between which places the spline is located:
            !
            ispg: DO isp = 1, spdim ! condition ca)
               IF (r < spxr(isp, 1)) CYCLE ispg ! distance is smaller than this spline range
               IF (r >= spxr(isp, 2)) CYCLE ispg ! distance is larger than this spline range
               ! at this point we have found the correct spline interval
               rz = r - spxr(isp, 1)
               IF (isp /= spdim) THEN
                  nsp = 3 ! condition ca
                  DO jsp = 0, nsp
                     erep = erep + scoeff(isp, jsp + 1)*rz**(jsp)
                  END DO
                  IF (dograd) THEN
                     DO jsp = 1, nsp
                        de_z = de_z + jsp*scoeff(isp, jsp + 1)*rz**(jsp - 1)
                     END DO
                  END IF
               ELSE
                  nsp = 5 ! condition cb
                  DO jsp = 0, nsp
                     IF (jsp <= 3) THEN
                        erep = erep + scoeff(isp, jsp + 1)*rz**(jsp)
                     ELSE
                        erep = erep + surr(jsp - 3)*rz**(jsp)
                     END IF
                  END DO
                  IF (dograd) THEN
                     DO jsp = 1, nsp
                        IF (jsp <= 3) THEN
                           de_z = de_z + jsp*scoeff(isp, jsp + 1)*rz**(jsp - 1)
                        ELSE
                           de_z = de_z + jsp*surr(jsp - 3)*rz**(jsp - 1)
                        END IF
                     END DO
                  END IF
               END IF
               EXIT ispg
            END DO ispg
         END IF
      END IF
      !
      IF (dograd) THEN
         IF (r > 1.e-12_dp) derep(1:3) = (de_z/r)*rv(1:3)
      END IF

   END SUBROUTINE urep_egr

END MODULE qs_dftb_utils

